Random death process for the regularization of subdiffusive anomalous equations 
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Subdiffusive fractional equations are not structurally stable with respect to spatial perturbations 
to the anomalous exponent (Phys. Rev. E 85, 031132 (2012)). The question arises of applicability 
of these fractional equations to model real world phenomena. To rectify this problem we propose 
the inclusion of the random death process into the random walk scheme from which we arrive at 
the modified fractional master equation. We analyze the asymptotic behavior of this equation, both 
analytically and by Monte Carlo simulation, and show that this equation is structurally stable against 
spatial variations of anomalous exponent. Additionally, in the continuous and long time limit we 
arrived at an unusual advection-diffusion equation, where advection and diffusion coefficients depend 
on both the death rate and anomalous exponent. We apply the regularized fractional master equation 
to the problem of morphogen gradient formation. 



Anomalous subdiffusion, where the mean squared dis- 
placement grows sub-linearly with time (a;^(t)) ~ 
where the anomalous exponent j/ < 1, is an observed 
natural phenomena 1]. It is seen in areas as varied as 
dispersive charge transport in semi-conductors ion 
movement in spiny dendrites protein transport on 
cell membranes [4]. In the classical paper Metzler, 
Barkai, and Klafter introduced the fractional Fokker- 
Planck equation (FFFE) that describes anomalous sub- 
diffusion of particles in an external field, F{x). This 
equation for the probability density p(a;, t) is written as 



a F{x) 

dx ksT 



-Qt - LFpp{x,t), where Lpp 
is the Fokker-Flanck operator, K^, is the anomalous diffu 
sion coefficient and 'DI~'^ is the Riemann-Liouville frac 
tional derivative of order 1 — v, defined as 
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p{x, T)dT 



(1) 



It was shown that the external field F{x) leads to a sta- 
tionary solution to the FFFE in the form of the Boltz- 
mann distribution [6]. However, in a recent paper [tI], we 
have demonstrated that this fundamental result is not 
structurally stable with respect to spatial variations of 
the anomalous exponent v. Small, non-homogeneous in 
space, variations of v destroy the stationary solution to 
the FFPE. In fact, even the simple one-dimensional frac- 
tional subdiffusion equation with constant anomalous ex- 
ponent and F{x) = 0, in the finite domain [0, L] with 
reflective boundary conditions, is structurally unstable. 
This equation should give a uniform stationary distribu- 
tion over the interval [0,1/] in the long-time limit. How- 
ever, if we use slightly non-uniform anomalous exponent 
v{x), the probability density p{x,t) will be completely 
different from the uniform distribution: as i — >■ oo it con- 
centrates at the point where ^{x) has a global minimum 
on [Ojij. We called this phenomenon anomalous aggre- 
gation [Sj] . Since it is impossible to have a completely ho- 
mogeneous environment, in which ly is uniform, the ques- 
tion arises as to whether the fractional equations with 
constant anomalous exponents are useful models for any 
real phenomena involving subdiffusion. This question is 



of great importance for the problem of a morphological 
patterning of embryonic cells, which is controlled by the 
distribution of signaling molecules known as morphogens 
To ensure robust pattern formation, the mor- 
phogen gradients must be structurally stable with respect 
to the spatial variations of environmental parameters in- 
cluding the anomalous exponent. Note that the unusual 
behavior of subdifFusive transport has been observed in 
an infinite system with two different values of anomalous 
exponents 12]. 

To rectify the structural instability involving unlimited 
growth of p{x, t), at the point of minimum of anomalous 
exponent v{x), we need a regularization of the fractional 
equations. The standard approach to regularize the frac- 
tional subdifFusive equations is to temper the power law 
waiting time distribution in a such way that the nor- 
mal diffusion behavior in the long-time limit is recovered 
(see, for example, [13] • In this case, by suppressing the 
power law behavior, the intrinsic characteristics of the 
anomalous process are lost. In this Letter we suggest a 
completely different approach, where we do not change 
the anomalous character, and retain the characteristics 
of the process. The main idea is to employ random death 
process and to 'kill' aging particles, for which the escape 
rate from the traps tends to zero as age tends to infinity. 

The main aim of this Letter is to show that as long a 
death process is introduced, together with a particle pro- 
duction, the stationary solution of the modified fractional 
master equation is structurally stable whatever the spa- 
tial variations of anomalous exponent might be. In par- 
ticular we use a regularized fractional master equation 
for the problem of morphogen gradient formation, which 
is a central topic of pattern formation in developmen- 
tal biology [13] ■ Here we deal with a discrete fractional 
master equation and its continuous approximation, cor- 
responding to a fractional Fokker-Planck equation. 

Let us consider a random walk of particles on a semi- 
infinite lattice with unit length. The particle performs a 
random walk as follows: it waits for a random time Tk 
at each point k before making a jump to the right with 
probability r{k) and left with the probability l{k). We 
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denote the residence time probability density function by 
ijj{k^T) = Pr{Tfc < r}, and assume it has the Pareto 
form 



'ip{k,T) 



iy{k)TQ 



a+i'(fc) ' 



where tq is a constant with the unit of time, and I'ik) is 
the spatially dependent anomalous exponent: < i^ik) < 
1. We assume that during the time interval {t, t + At) at 
point k the particle has a chance 9{k)At + o{At) of dying, 
where d{k) is the death rate {0{k) > 0). 

We denote by p{k, t) the average number of particles at 
the point k at time t. The anomalous subdiffusive master 
equation with the death process can be written as 

+ b{k + i)e-''('=+i)*pi-''(^-+i) [p{k + 1, t)e^('=+i)*] 
-e{k)p{k,t), k>2 (2) 



where 



aik) = 



r(i-Kfc)K 



h{k) = 



m 



r(i-Kfc)K 



u(k) 



are the anomalous rate functions. This fractional equa- 
tion can be derived from a number of standpoints (see, 
for example, [13] )• In this equation the anomalous ex- 
ponent depends on the state, which is crucial for what 
follows. For the case of constant anomalous exponent v, 
this reaction-transport equation and its continuous ap- 
proximations were considered in 

To ensure the existence of stationary structure in the 
long time limit, we introduce the constant source term g 
at the boundary of the semi-infinite lattice (fc = 1). This 
is crucial for the problem of morphogen gradient forma- 
tion, where g models a localized source of morphogens 



We assume that the boundary is reflective, so we 



have the equation for p{\,t) 



dp{l,i) 
dt 



b{2)e 



-e(2)t^i-f(2) 



vr''>[p{2,t). 



a{l)e- 



=9(1)41 



ep{l,t)+g. (3) 



Note that any nonlinear function g{p) can be included in 

Without the reaction { 9 = Q) the fractional mas- 
ter equation ^ with constant anomalous exponent v is 
structurally unstable in the long time limit. The station- 
ary solution pst{k) ~ limt^ oo p{k,t) can be found from 
as 



fc-2 



Pst 



b{j + iy 



(4) 



where 



oo k—1 



k=2 j = l ' 



(5) 



provided the sum is convergent. However, when the 
anomalous exponent is not constant, the asymptotic be- 
havior is completely different. Consider the point M, at 
which the anomalous exponent is at a minimum v{M) < 
v(k), Vfc ^ M. Then, one can show that 



p(M,i)^l, p{k,t)^Q, 



t — > oo. 



(6) 



As stated earlier, the main aim of this Letter is to reg- 
ularize the fractional Master equation with the addition 
of the random death process. To this end, it is convenient 
to rewrite the fractional master equation as 



dp{k,t) 
dt 



= -I{k, t) + I{k - 1, - 0{k)p{k, t), 



k>2 
(7) 



where /(fc, t) is the total flux of cells from k to k + 1 

I{k, t) = a(fc)e-«('=)*I?^''^"^ [p(fc, i)e''('=)*] 

^ b{k + i)e-«('=+i)*2?i-^(''^+i) [p{k + 1, t)e''(''^+i"] 



(8) 



The flux I{k,t), in Laplace space takes the form 

I{k, s) = a(fc)(s + 0(fc))i-"('=)j5(fc, s) 

- b{k + l)(s + 6{k + l))i-''('=+i)j5(A: + 1, s). 



(9) 



From here we can find the stationary flux Ist{k) = 
lims_).o sl{k, s) as follows 

ht{k) = a*{k)p,t{k) ~b*{k + l)pst{k + 1), 

where 

a*{k) = a{k) [eik]]^'"^''^ , b*{k) = b{k) [0{k)]^-''^''^ 

and Pst{k) = lims^o sp{k, s). The main feature of this 
stationary flux is that it has Markovian form; but the 
rate functions a*{k) and b*{k) depend on the anomalous 
rate a(fc), b{k), random death rate 9{k), and anomalous 
exponent This unusual form of stationary flux is 

because of the non-Markovian character of subdiffusion. 

Let us flnd the stationary distribution Pst{k) for the 
simple case where 9 is constant. In the long time limit, 
at the boundary, we then have the following condition: 



htii) = g-Opstii). 



(10) 



Similarly, we have the condition at the location k = 2 
Ist{'2) — -fst(l) — 0Pst{'^)- We are able to obtain a general 
expression for the stationary flux at location k 



Ist{k) 



k 



Pstij) 



(11) 



3 



This has a very simple physical meaning: that as t — cxd, 
Istik) tends to the difference between the production rate 
and the sum of death rates at all states from the bound- 
ary up to k. It is clear that as — >■ oo, the stationary 
flux Ist{k) 0, since in the stationary state g should be 
equal to total death rate 



.9 = d^Pst{j)- 

We obtain 

b{k + l)9-'^''+^^p,t{k + 1) = a{k)e-'^''^p,t{k) 



(12) 
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This equation allows us to find Pst{k) for all k. For the 
symmetrical random walk for which a{k) = b{k) = a and 
v = const, we have 



Pstik + 1) = - - [ I - J2Mj) I • (13) 



Now let us obtain the fractional Fokker-Planck equation 
with the death process, as the continuous limit of the 
master equation We change the variables A: — >■ a;, 
k zL 1 X zL I and take the limit / — > to obtain 



dp 



_d_ 

dx 

dx'^ 



I {a{x) - b{x)) e-«(-)*pi-''(^) [p(x, t)e'^^^'] 

^(a(:r)+6(x))e-«(-)*Pr''(^)[p(a:,t)e«(-)*] 
^9ix)p{x,t). (14) 



From this equation we obtain for Pst{x) — \\m.t^oo p{x, t) 
the stationary advection-diffusion equation 

^ [vl {x)pst{x)] + ^ [Dl {x)p,t{x)] = 0{x)p^t{x), 



dx 



where vf, {x) is the drift, and (x) is the generalized 
diffusion coefficient defined as 



I {r{x) - Ijx)) [eix)]^-"'^''^ 
r(l - iy{x)){Tor(.-) - 



Dl {x) = 



p [e{x)] 



2r(l-z.(x))(ro )"(-)■ 



This result means that in the long time limit, subdifFusion 
with the death process becomes standard diffusion with 
nonstandard drift {x) and diffusion coefficient Df,{x). 
Both of them depend on the death rate d{x) and the 
anomalous exponent iy{x). This is due to non-Markovian 
character of subdiffusion. It has been found in lldl that 



the non-Markovian behavior of subdiffusion leads to an 
effective nonlinear diffusion. Note that the drift term 

(x) plays an essential role in chemotaxis, since vf, (x) ^ 
where C is the chemotactic substance. Therefore 
the dependence of chemotactic term of the degradation 
rate 9 can be of great importance for the problem of cell 
aggregation [i,[ia|20]. 

Let us consider a random walk with a constant drift 
vf, — ~v, diffusion D^, and degradation rate 9. Then 



dPst{x) gd'^pstix) 



dx ^ dx^ 

The solution is the exponential profile 



9pst{x)^Q. (15) 



Pst{x) = Aexp 



+ ^y-^ + i.Dl9 



(16) 



where A can be found from the condition g = 
^ Pstix)dx : 

g(v + ^v^ + 4Df,9) 

When — 0, we have a morphogen profile obtained in 
111: 



Pst{x) 



exp 




(18) 



We now simulate the fractional master equation with ran- 
dom death process, (O, using Monte Carlo techniques. 
Throughout this we let tq = 1, so that this is the unit 
of time for the simulation; we take g = 1, so that we 
have a constant birth rate of one particle per unit time. 
The first particle begins a random walk at = 1, such 
that at each point k waiting times are distributed as 



ijj{k,T) 



and jump probabilities to the left 



(to+t)1 + -<'=) ' 

and right from each point k are r(fc) and l{k) respectively. 
A particle completes a random walk from when it is pro- 
duced until the terminal time t = T, or until its random 
time of death exponentially distributed as ijjDit) = 9e~^* . 
This death rate is equivalent to a spatially invariant, con- 
stant death rate in . Also note that unlike the wait- 
ing time, the death time is not renewed when the particle 
makes a jump. The practical issue of having particles be- 
ing produced and dying is dealt with in the following way. 
The first particle in the simulation begins at time t = Q, 
and completes its random walk as described above; the 
second particle begins at t = 1, because tq = 1, and 
completes its random walk; and so on until time t = T . 

Firstly let us consider the symmetrical random walk, 
where r[k) = l{k) = i, v{k) = 0.5, and 9 = lO^^. The 
figure FIG.[l]shows the corresponding stationary density 
made up from 10^ realizations of the random walk at 
time T = 10^. 
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FIG. 1. Stationary profile for tiie symmetric fractional master 
equation where r{k) = l{k) = |, v{k) = const. = 0.5, to = 1, 
and e = 10"^ 
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FIG. 2. Stationary profile for the symmetric fractional master 
equation, with a pertubation to the anomalous exponent at 
k = 8. u{k 7^ 8) = 0.5, 1^(8) = 0.4 



We can see that our simulation is in agreement with the 
analytical values calculated from the recurrence relation 

Next, we show that the model is robust to non- 
homogenous spatial perturbations in the anomalous ex- 
ponent. Analogously to the simulation we presented in 
the previous work [7[ , we introduce a small perturbation 
to the anomalous exponent at one point in the space: all 
states have v = 0.5 except for fc = 8, which has v = 0.4. 
From FIG. [2] we can see that although we observe a 
change to the stationary profile around point k = S, the 
stationary profile is structurally stable and exponential 
in character. We stress the importance of the death pro- 
cess in regulating the behavior of the process to ensure 
stability. Whereas in our previous work, we showed that 
even a small perturbation in the anomalous exponent like 
this would lead to a breakdown in the stationary density. 
Additionally, we considered a non-symmetrical random 
walk which leads to a drift and found that the profile is 
stable. 

In summary, we have suggested a new regularization 
of the subdiffusive fractional master equation by using 
the random death process. The fundamental feature of 
this approach is that unlike the previous regularization, 
which loses the anomalous characteristics, we are able to 
retain dependence on the anomalous exponent. We find 
the stationary flux of the particles has a Markovian form. 



with unusual rate function depending on the anomalous 
rate functions, the death rate, and the anomalous expo- 
nent. We have shown that the long-time and continuous 
limit of this regularized fractional equation is the stan- 
dard advection-diffusion equation that, importantly, is 
structurally stable with respect to spatial variations of 
anomalous exponent v. We have found that the effective 
advection and diffusion coefficients, and Df,, are in- 
creasing functions of the death rate 9: vf, ^ ^ 6^^" . 
We have applied a regularized fractional master equa- 
tion and modified fractional Fokker-Planck equation to 
the problem of the morphogen gradient formation. We 
have shown the robustness of the stationary morphogen 
distribution against spatial fluctuations of anomalous ex- 
ponent. 
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